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Abstract 

We present a multigrid scheme for the solution of finite-element 
Hartree-Fock equations for diatomic molecules. It is shown to be fast 
and accurate, the time effort depending linearly on the number of 
variables. Results are given for the molecules LiH, BH, N2 and for the 
Be atom in our molecular grid which agrees very well with accurate 
values from an atomic code. Highest accuracies were obtained by 
applying an extrapolation scheme; we compare with other numerical 
methods. For N2 we get an accuracy below 1 nHartree. 
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1 Introduction 



For benchmark calculations in Hartree-Fock it is necessary to have a very 
good accuracy. This may be achieved through numerical methods like finite 
differences |21EllllEllSllZlElinilIIIllIIllI2llinilIl!. Recently, the accuracy 
has been pushed to the sub-//Hartree range by employing more than 10000 
points in a 2-dimensional finite difference grid^l4j. Alternatively a finite- 
element program was developed which gave quite accurate results [TTl HH] 
with less points. However, the time effort of this method scales quadratically 
in the number of unknowns. 

In recent papers by Kopylow et al. |2H |22] it was shown that a multi- 
grid approach for solving finite element equations applied to the Kohn-Sham 
equations of density functional theory gives fast and accurate results. 

Extrapolation methods are a tool to gain better accuracy from 

a sequence of values for a given property for which the asymptotic behaviour 
is known. They have been successfully applied to MP2 correlation ener- 
gies for closed shell atoms j2Sl where several orders of magnitude have been 
gained. More recent is the application to one-electron Dirac-FEM solutions 
for diatomics [2Ij which yields also several orders gain in accuracy. 

In this paper we present a scheme which utilizes multigrid techniques 
|19L OU] and scales linear in the number of variables. Extrapolation methods 
are employed in addition in order to further improve the accuracy. 

In part |21 we describe the Hartree-Fock approximation used. In part El 
our coordinate transformation is given. In part 0] we state the discretization 
of the Hartree-Fock equations by the finite-element method. And in part El 
we explain the multigrid scheme that is used. In part (HI we give results both 
directly calculated and extrapolated for Be, BH, LiH and N2 and compare 
with other work[21 ITllTij. 
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2 Hartree-Fock Method 



Microscopic physical systems like molecules are described by wavefunctions 
whose behaviour is governed by the Hamiltonian of the system. Often one 
is only interested in the energy levels of the system which are given through 
the eigenvalue equation 

= E<iJ (1) 

The eigenfunction \& is a many-body wavefunction. Equation^ can be solved 
analytically for a few physical systems only. For others it is necessary to 
make approximations like the Hartree-Fock method where a single determi- 
nant serves as the ansatz function for the many-body wave function(Slater 
determinant). For electronic systems the Hamiltonian is 

N N 

H = Y.T^ + T.V^ + T.Vee{\n-r,\) (2) 
i=l i=l i<j 

where Tj denotes the kinetic energy operator — of a single electron, Vi = 
Vexti^i) is a given external potential and Vee the interaction potential between 
electrons. The variation of the single particle wave functions ^^(orbitals) in 
the Slater determinant results in the Hartree-Fock equations: 

{f + V;,i(f) + VD^r{r)} 0,(r) - E yd^<l>^i^ = ^jM^ (3) 

i 

AV,.{f) = -4np,,{f) (4) 

AVD^ri^ = -47rp(f) (5) 

where p{f) = p{f, r) is the density, the diagonal part of the density matrix 
p{f,r') = Z^i 0i and Pji{r) = (f)j{f)(j)*{f) the exchange density. We 

now define gl = J2i Vji{r)(t>i{r) and Hiok = f + Vext{r) + Voirir). EquationEl 
then takes the form 

Hiok(t)j{r) -gi = ej(j)j{r) (6) 

Equations 131101 are solved iteratively (SCF iteration) in the following way. 
First we take a trial set of functions compute Voir and Vji and get Hiok 
and gi. Then one calculates an (approximate) eigenvalue if 

e,- = j (f)*{Hiok(pj - gi)d^r (7) 
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and solves directly the inhomogeneous equation system 

{Hiok-~e,)<Pj = gi,J = l,...,N (8) 

in order to get new wave functions {(pi}. Then these new wave functions are 
orthogonalized from the functions with the lowest to the highest eigenvalues 
ii of occupied states(Gram Schmidt method). For these wavefunctions the 
total energy is computed as the expectation value for the corresponding Slater 
determinant. From these new Voir, Vji and Hiok, gi are computed and 
equations [7[|H1 solved again. This process is re-iterated till a wanted accuracy 
in the total energy or the energy levels of the orbitals is obtained. In order 
to get stable and fast convergence the direct potential is mixed with the old 
one: Voir = Pmix * V^[i + (1 - Pmix) * Ygf^"; pmix = 0.95. 
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3 Diatomic molecules 



Diatomic molecules have axial symmetry, commutes with H and thus 
leads to good quantum numbers nij = Izj- We take a restricted Hartree-Fock 
approach and demand the symmetry of the single-particle wavefunction to 
have axial symmetry also. This leads to the ansatz in cylindrical coordinates 
for the orbitals. 

0,(r,^,^)=/,(r,^)e^"^^-^ (9) 

The two-center point nucleus Coulomb singularities are best described by 
elliptic hyperbolic coordinates: 

« = ^ (10) 

V - (11) 

where Vk is the distance to th k-th. nucleus(k=l,2) and R the internuclear 
distance. These coordinates remove the singularity of the Coulomb poten- 
tial which is necessary for a high order convergence behaviour of the finite- 
element method. For the computation of best energies in closed shell systems 
the prolate spheroidal coordinates are favorably used which emerge after a 
singular coordinate transformation (the back transform is singular at the nu- 
clear centers, i.e. at ^=l(s=0), -q — rkl{t — 0, tt)). 

i = coshs (12) 
rj — cost (13) 

This singular transform improves the analytic properties of polynomial ansatz 

functions in s.t considerably in a finite-element approach even though it 
leads to a problematic behaviour at the inner boundaries ^=1 (s=0) and 
?7=±l(t=0,7r) respectively. In FEM this difficulty is adequately handled by 
an open boundary treatment. On the outer boundary we use a closed bound- 
ary with boundary values 0. 
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4 Finite element method 



In the finite element method a variational formulation is generally the start- 
ing point. In our case, the variational integral corresponding to equation |S1 
is: 

/ = J cl>;{^{Hiok - ~ej)<Pj - 9i)d'r (14) 

In the finite-element method(FEM) the space is subdivided into several 
subspaces called elements on which locally defined formfunctions Ni are used 
as approximation. Via the FEM ansatz 0(f^ = J2i CiNi{r) and variation with 
respect to the Cj the following equation results: 



lok,kl 



Si) 



yx,k 



lok,kl 

yx,k 



Ski 



(15) 

(16) 

(17) 

(18) 
(19) 

At the outer boundary the wavefunctions are set to zero(closed bound- 
ary), at the inner boundaries (symmetry axis) the wavefunctions may take 
any values(open boundary). Analogously there exists a variational integral 
for equation HJ 



- ejSki)c\- 

Nl{r)Hi,uNi{r)dh 
Nl{T)gidh 
Nl{f)Ni{r)d^r 



I 



(20) 



A corresponding integral exists for the direct potential from equationElwhich 
then can be treated in the same way. This leads after insertion of the FEM 
ansatz to the linear inhomogeneous equation systems 



with Dik = j VN;{r)VNk{r)d^r 
and pi^ = An[ pij{r)Nl{r)d^r 



(21) 

(22) 
(23) 
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The outer boundary values for the potentials Vij are computed from a 
multipole expansion of the densities pij up to the fifth order. The ansatz 
functions Nk are complete polynomials of order p (p < 7) in s,t and thus 
complicated transcendental functions of the real space vector r; p is called the 
order of the elements and is normally taken to be 4. At the element bound- 
aries continuity at the grid points is demanded which due to the serendip 
properties of polynomials leads to continuity over the whole element bound- 
aries. It can be shown that for large numbers of points n the error due to 
the finite number of points is proportional to ^pQ. 
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5 Multigrid method 



In order to solve the FEM matrix equations we use a multigrid scheme. 
The multigrid method ^3 120] combines the good smoothing behaviour of 
iterative methods(e.g. Gauss-Seidel, CG) with an effective elimination of 
long-range errors which are treated at coarser grids. In diagram ^ we show 
the MG scheme we use. First we solve directly on the coarsest grid(E). Then 
this solution is prolongated to the next finer grid(P5). Prolongation is an 
interpolation step to the finer grid values where we use all ansatz functions 
over each coarse element. If we have a solution from an older scf cycle we 
compare the defects of both and take the one with the smaller one(V). If 
this vector is not converged we do a V-cycle: First we restrict the defect to 
the next smaller gTid{RD)- The restriction step is done by the transpose of 
the interpolation matrix Ps,Rd = -Pj and is thus of the same (high) order 
as the interpolation. If this is not the coarsest grid, we smooth and restrict 
the resulting defect to the next coarser grid(S,i?z))- This is repeated till 
the coarsest grid is reached where we solve directly for the defect (D). The 
solution of the defect is prolongated to the next finer grid and added to the old 
approximation vector(Pc';+)- This new vector is then smoothed(S). This is 
repeated till the level from which the V-cycle started is reached. If the vector 
is not converged well enough we repeat the V-cycles till sufficient convergence. 
Then we prolongate the solution vector(P5) to the next finer grid and repeat 
the whole procedure. In our computations we used 15 conjugate gradient 
steps in each smoothing. On the average we needed 2 V-cycles on the finest 
grid per scf iteration for every state and potential Vij. 
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Figure 1: multigrid algorithm 
F: finest grid; C: coarsest grid 




repeat until convergence 
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6 Results and Discussion 



We first present values for BH, LiH and Be using different grids ranging from 
545 to 16641 points. Tliese results were extrapolated using both inverse 
power and rational geometric extrapolation which are described in Flores 
and Kolbj2ni- For large numbers of points n the energy obeys the following 
formula: 

E{n) = E^ + - (24) 
nP 

where E(n) is the energy for a given number of points, Eoo the exact value, 
p the polynomial order and C a constant. The same relationship holds for 
all properties except the densities at the nuclei. For the extrapolation of the 
densities at the nuclei we took p = 3 for the leading term of the asymptotic 
expansion. Fitting of three subsequent densities to the formula D{n) = 
Doo + ^ gave a convergence of the parameter p close to 3. 

In table ^ the properties of BH(R=2.336) are given for different num- 
bers of points and the extrapolated values together with the finite difference 
results of Laaksonen et al.|2] and, for the total energy and electron levels, 
Kobus[7j. Compared with those of Laaksonen et al. our directly computed 
results are generally more accurate up to three digits. The only exception 
are the densities at the nuclei. For the hydrogen nucleus we are as accurate 
as Laaksonen et al. and at the heavier nuclei our directly calculated results 
are about two digits less accurate, however by extrapolation we gain 4 dig- 
its. This reflects the fact that the finite element method optimizes integral 
properties not point-like ones. In most cases our values are different from 
Laaksonen et al.'s though mostly only in the last digit given by them. 

Our results for 16641 points are as accurate as those of Kobus, who had 
more than three times the number of points(52441). And our extrapolated 
results should be clearly more accurate. 

For LiH(R=3.015) the extrapolated values are given in table El Our 
results agree with the more recent numbers of Kobus but give 2-3 more 
digits. Be is computed with the two-center grid where the Be atom is placed 
in one center, the other center, R apart, is empty(dummy center). We can 
compare with the results of the atomic GRASP code[28.;_. The total energy 
is given for 390 points, the energy levels, due to convergence problems, for 
220 points. For the total energy our result disagrees slightly from the Kobus 
values but agrees with the GRASP values which shows our values to be more 
accurate. For the energy levels the GRASP results agree with both but has 
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to few digits to differentiate between them. 

In figure |21 the computation times for the solution routines for different 
numbers of grid points are given. It shows a linear dependence on the number 
of points. The same is true for the total time per scf iteration. This allows 
us to use rather large numbers of grid points. 

In table 121 the convergence of the total energy of N2(R = 2.068) with 
respect to the number of points is given up to 148225 points. Here we have 
the expected convergence behaviour of the finite element method where the 
leading error term for a high number of grid points is proportional to ;^(p 
being the order of the polynomials used). The extrapolated values were 
computed for one sequence without and one with the 148225 point value. 
Both are identical apart from the self-consistency error. In order to test the 
accuracy of the results we took the truncation parameter d=18 a.u. instead of 
25 a.u. in order to test whether the results were dependent on the properties 
on the outer boundary, d is the distance between the point on the outermost 
ellipse ^max = const, and a focal point, if the distance is taken perpendicular 
to the symmetry axis. The values show a slightly faster convergence because 
of the higher density of points in the inner region but converge to the same 
result for the higher number of points. 

In table |3] the results for 6th-order polynomials with d = 40 a.u. and 
7th-order polynomials with d = 18 a.u. are shown. No effect from bound- 
ary values or from the different orders can be seen up to InHartree. The 
difference of the various orders is probably due to the bigger truncation er- 
ror accumulation for the higher orders. This shows our results for the total 
energy of N2 to be accurate to at least 1 nHartree, and presumably up to 
two more digits. The finite difference result of Kobus et al. Jl] has an error 
of 13 nHartree for his 793 x 793 grid lying well in the sub-yuHartree level of 
accuracy as was claimed. This accuracy can be reached with our standard 
method with only 37249 points. For higher orders p it takes only 9409 (p=6) 
or 3613(p=7) points. In comparison to the finite difference scheme of Kobus 
et al.[14j we can achieve an accuracy which is by al least 2 orders of magni- 
tude better with much smaller numbers of grid points and thus remarkably 
small computational times for our high precision benchmark results. The 
same holds true with respect to the old finite-element results of Heinemann 
et al. [IH], where we gain 5-7 digits(see table E)). It should be noted that 
our computation with the highest number of points took less than a day 
on an ordinary personal computer. At last, we want to point out that ex- 
trapolation schemes have to be applied judiciously. In figure |2] the relative 



11 



errors of the total energies are given for different grids. Unlike the expected 

convergence oc ;^(p=4) one gets an alternating order parameter p. Closer 
inspection showed that the grids have alternatingly a different geometry. If 
grids with the same geometry are taken one gets p — > 4 and correspondingly 
good extrapolation values. 
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Table 1: Results for BH(R=2.336): total energy, energy levels 
moments of order L and densities at the nuclei; all in a.u. 



multipole 



points 


total energy 


ei 


62 


es 


625 


-25.13098274805960 


-7.68616024068801 


-0.648192901104539 


-0.34 8418481279812 


1089 


-25.13142119823359 


-7.68622122038152 


-0.648187699192955 


-0.348423192223713 


2401 


-25.13159271731769 


-7.68626555785907 


-0.648187293511312 


-0.348423751988741 


4225 


-25.13159798871873 


-7.68626718236676 


-0.648187268672215 


-0.348423777833311 


9409 


-25.13159867258947 


-7.68626736927568 


-0.648187265330006 


-0.348423781443339 


16641 


-25.13159869928299 


-7.68626737702293 


-0.648187265204238 


-0.348423781582957 


extrapol 


-25.13159870231 


-7.68626737794 


-0.6481872651901 


-0.3484237815982 


Kobus^ 


-25.13159870 


-7.686267370 


-0.648187256 


-0.348423779 


Laak.|2] 


-25.131647 


-7.686283 


-0.648190 


-0.348426 



points 


L=l 


2 


3 


4 


625 


-5.3526511173 


12.1863258126 


-15.6409316245 


24.8481986257 


1089 


-5.35 24735227 


12.1862527191 


-15.6411035350 


24.8492335950 


2401 


-5.35 2 4679377 


12.186 2405731 


-15.6410934947 


24.8492589014 


4225 


-5.3524669842 


12.1862400311 


-15.6410933554 


24.8492609410 


9409 


-5.3524669579 


12.1862399666 


-15.6410933548 


24.8492612346 


16641 


-5.3524669569 


12.1862399640 


-15.6410933546 


24.8492612446 


extrapol 


-5.3524669568 


12.1862399637 


-15.6410933546 


24.8492612457 


Laak.0 


-5.352466 


12.18621 


-15.64103 


24.84888 



points 




pirn) 


625 


69.179642371 


0.46683681664 


1089 


70.644543782 


0.46741415376 


2401 


71.536579057 


0.46 754727884 


4225 


71.660901149 


0.46 755872297 


9409 


71.691078340 


0.46756104352 


16641 


71.693793197 


0.46756122952 


extrapol 


71.694413 


0.4675612700 


Laak. [2^ 


71.69451 


0.467561 
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Table 2: results for LiH,Bc: total energy, energy levels e^, multipole moments 
of order L and densities at the nuclei; all in a.u. 



LiH 
R=3.015 


extrapolated values 


Laaksonen et al. 


Kobus 


energy 


-7.987352237228 


-7.987354 


-7.987352237 


ei 


-2.445233713306 


-2.445234 


-2.4452337133 


e2 


-0.301738270249 


-0.301738 


-0.3017382702 


L= 1 


-0.65318943587 


-0.653190 




2 


7.12821973712 


7.128219 




3 


-2.90955527116 


-2.909556 




4 


16.275742582 


16.02756 




p{rLi) 


13.789722803 


13.789729 






0.37406093101 


0.374061 





Be 


extrapolated values 


Laaksonen et al. 


Kobus 


GRASP 




R=2.00 




R=2.00 




energy 


-14.573023168305 


-14.5730226 


-14.573023170 


-14.573023168 




-4.732669897448 


-4.7326689 


-4.732669898 


-4.7326699 


62 


-0.3092695515724 


-0.30926957 


-0.3092695522 


0.30926955 



Table 3: Total energy of N2(R=2.068) 



points 


total energy(d=25) 


total energy(d=18) 


625 


-108.988278969512 


-108.990118887918 


1089 


-108.992300729719 5 


-108.992661731070 


2401 


-108.993744526847 


-108.993775395838 


4225 


-108.99381783520 


-108.99382035051 


9409 


-108.993825276408 


-108.993825404018 


16641 


-108.993825597923 5 


-108.993825610993 


37249 


-108.99382563334 


-108.99382563387 


66049 


-108.99382563467 


-108.99382563472 


extrapol 


-108.99382563482 


-108.99382563482 


148225 


-108.99382563481 


-108.99382563482 


extrapol 


-108.99382563482 


-108.99382563482 


Kobus et al. 


-108.993825622 
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Figure 2: time of solution routines for 1 scf iteration 

N2 

solid: for potentials; dashed: for wavef unctions 



50 



5000 10000 

number of grid points 



15000 



20000 



Tabic 4: Total energy of N2 for different orders and d 



points 


total energy (p=6,d=40) 


points 


total energy (p=7,d=18)) 


1369 


-108.9938095840 


1849 


-108.9938248739 


2401 


-108.9938243325 


3613 


-108.9938256311 


5329 


-108.9938256092 


7225 


-108.9938256350 


9409 


-108.9938256343 






21025 


-108.99382563487 






37249 


-108.993825634866 
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Table 5: results for N2: total energy, energy levels ej, multipole moments of 
order L and densities at the nuclei; all in a.u. 



points 


extrapolated values 


Heinemann et al. [T^ 


energy 


-108.99382563482 


-108.993826 




-15.68186695242 


-15.681867 




-15.67825164397 


-15.678252 


es 


-1.473422499578 


-1.473423 


£4 


-0.778076815628 


-0.778077 


£5 


-0.6347931345534 


-0.634793 


ee 


-0.6156250666967 


-0.615625 


L= 2 


15.908084537079 




4 


23.3942874333 






205.3983861 
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Figure 3: convergence for grids with different geometry 



black and white circles denote different geometries 




1000 10000 100000 

number of gridpoints 



17 



7 Acknowledgement 

One of us(O.B.) acknowledges financial support of the Deutsche Forschungs- 
gemeinschaft (DFG) . 



18 



References 

[1] Strang G, Fix G 1973 An Analysis of the Finite Element Method (En- 
glewood Cliffs: Prentice-Hall) 

[2] Laaksonen L, Pyykko P, Sundholm D 1983 Chem. Phys. Lett. 96 1 

[3] Pyykko P, Sundholm D, Laaksonen L 1987 Mol. Phys. 60 597 

[4] Sundholm D 1988 Chem. Phys. Lett. 149 251 

[5] Pyykko P 1989 Numerical Determination of the electronic Structure 
of Atoms, Diatomic and Polyatomic Molecules ed M Defranceschi, J 
Delhalle(Dordrecht: Kluwer) p 162 

[6] Pyykko P, Sundholm D, Laaksonen L, Olsen J 1991 Europhysics Con- 
ference on Computational Physics ed A Tenner (Singapore: World Sci- 
entific) p 455 

[7] Kobus J 1993 Chem. Phys. Lett. 202 7 

[8] Kobus J 1994 Comp. Phys. Comm. 78 247 

[9] Moncrieff D, Kobus J, Wilson S 1995 J. Phys. B 28 4555 

[10] Kobus J, Laaksonen L, Sundholm D 1996 Comp. Phys. Comm. 98 346 

[11] Moncrieff D, Kobus J, Wilson S 1998 Mol. Phys. 93 713 

[12] Kobus J, Moncrieff D, Wilson S 1999 Mol. Phys 96 1559 

[13] Kobus J, Moncrieff D, Wilson S 2000 Mol. Phys 98 401 

[14] Kobus J, Quiney H M, Wilson S 2001 J. Phys. B 34 2045 

[15] Kobus J, Moncrieff D, Wilson S 2001 J. Phys. B 34 5127 

[16] Kobus J, Moncrieff D, Wilson S 2002 Mol. Phys 100 499 

[17] Heinemann D, Pricke B, Kolb D 1988 Phys. Rev. A38 4994 

[18] Heinemann D, Rosen A, Fricke B 1990 Physica Scripta 42 692 



19 



[19] Hackbusch W 1985 Multigrid Methods and Applications (Berlin: 
Springer) 

[20] Brandt A 1984 Multigrid Methods: 1984 guide with applications to fluid 
dynamics (Bonn: Ges. f. Mathematik u. Datenverarbeitung) 

[21] V Kopylow A, Heinemann D, Kolb D 1998 J. Phys. B 31 4743 

[22] V Kopylow A, Kolb D 1998 Chem. Phys. Lett. 295 439 

[23] Brezinski C, Zagha M R 1991 Extrapolation Methods. Theory and Prac- 
tice (Amsterdam: North-Holland) 

[24] Walz G 1996 Asymptotics and Extrapolation (Berhn Mathematical Re- 
search) 88 (Berlin: Akademie) 

[25] Stocr J, Bulirsch R 1992 Introduction to Numerical Analysis 2-nd Edi- 
tion (Berlin: Springer) 

[26] Flores J R, Kolb D 1999 J. Phys. B 32 779 

[27] KuUie O, Kolb D 2001 Eur. Phys. J. D 17 167 

[28] Surzhykov A, private communication 



20 



